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SUMMARY 

Higher-order parametric discrete elements are a significant modeling 
advance over similar elements with straight- sided triangular or quadri- 
lateral planforms. However, the addition of such discrete elements to 
NASTRAN poses significant interface problems with the Level 15. 1 assembly 
modules and geometry modules. The present paper systematically reviews 
potential problems in designing new modules for higher-order parametric 
discrete elements in both areas. An assembly procedure is suggested that 
separates grid point degrees of freedom on the basis of admissibility. New 
geometric input data are described that facilitate the definition of surfaces 
in parametric space. 


SYMBOLS 

k 

C Denotes continuity through k derivatives 

ff a p The partial derivative of f with respect to a and /3 

df Total differential of f 

R Cylindrical radius 

S Arc length 

u a Elastic displacement in the curvilinear coordinate direction a. 

3 

u Elastic displacement normal to the midsurface 

UA Grid point displacement parameters required for admissibility 

*This work was performed under the sponsorship of the McDonnell Douglas 
Astronautics Company Independent Research and Development Program. 
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UH 


Grid point displacement parameters not required for admissibility 
x 1 Grid point coordinates in a reference coordinate system 

Q l Rotations about the x 1 coordinate directions 

(|, v) Patch parameters analogous to curvilinear coordinates 

Subscripts: 

FP Flat plate 

CP Cylindrical panel 


INTRODUCTION 

In the main, joining problems arise with higher-order discrete elements 
because the additional grid point degrees of freedom contain terms directly 
proportional to element strains. A complete one-to-one joining of displace- 
ment parameters between geometrically similar elements implies a strain 
continuity that is erroneous if the elements are of different materials. 

A complete one-to-one joining of displacement parameters between geo- 
metrically dissimilar elements implies a strain discontinuity that is erron- 
eous if the elements are of similar materials. An example of this latter 
behavior will be given in which the in-plane displacement gradients are 
erroneously linked between cylindrical panel elements and flat plate elements. 
To obtain the correct solution these parameters must either be allowed to 
vary independently or be joined by a constraint equation for in-plane strain 
continuity. In general, any variational problem can be solved by using only 
the minimum constraints required to produce an admissible displacement 
field for the Ritz procedure. As a practical matter this approach is of little 
help if used uncritically since it can add many unnecessary degrees of free- 
dom. At the intersection of four MDAC parametric discrete elements, for 
example, there are 48 degrees of freedom that are reduced by admissibility 
constraints to 22 independent degrees of freedom. When strain continuity 
constraints are applied these are reduced to 12 independent degrees of free- 
dom. In this case there is nearly a 50 percent reduction when strain contin- 
uity is valid. Any modification of the NAS TRAN assembly modules to process 
higher -order discrete elements must be flexible enough to take advantage of 
this situation if it is to be efficient. 
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The basic geometric entity used by NAS TRAN is the grid point. The 
basic geometric entity need for parametric discrete elements is a mapping of 
two surface coordinates into points on the surface in three dimensions. This 
mapping, called a patch, is approximated locally by interpolation functions. 
These functions must be input to NAS TRAN in order to generate element 
matrices (stiffness, etc. ). It is, of course, feasible to input the patch for 
each element directly as part of the property data for the element. This has 
the obvious disadvantage of requiring a great deal of input data; up to 48 items 
for a bicubic patch. To reduce input data requirements, the MDAC para- 
metric plate element program uses the boundary curve for the entire plate to 
generate patches for each discrete element once a topological mesh has been 
specified (Reference 1). The present paper considers new geometric input 
data for NAS TRAN to facilitate the introduction of parametric discrete elements, 

ADMISSIBLE DISPLACEMENT FIELDS 

Admissibility conditions for discrete element displacement functions are 
an especially important topic for higher -order discrete elements. In class- 
ical variational mechanics the material properties are usually assumed either 
constant or continuously differentiable functions of the spatial coordinates. 

This leads to simple smoothness requirements based on the order of the dif- 
ferential operator in the equilibrium equations (Reference 2). The displace- 
ment u^ in a homogeneous plate bending problem, for example, must be 
in the interior and on a free edge. The latter condition is a consequence 
of the natural boundary condition for shear. These are of course conditions 
on the continuum displacement solution that the discrete element model must 
converge to in the limit. Convergence is measured by an energy norm and 
it is the existence of this norm that sets the admissibility conditions for the 
piecewise polynominal displacement function formed by assembling individual 
discrete elements. Returning to the plate example, the energy norm is 
derived from the strain energy density which involves at most second deriva- 
tives of u^. As long as the discrete element displacement field is at least 
C * between elements the energy norm is well defined. Physically this con- 
dition corresponds to the absence of a hinge between plate elements and it is 
imposed at the grid points to assemble or build a discrete element model of 
a plate structure. How closely the assumed displacement functions approach 
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this condition between nodes is a problem in approximation theory that is 
intimately related to the question of completeness. This is another issue 
entirely and for the present discussion completeness will be assumed. Again 
returning to the plate example if there are no line moments between elements 
and the material is continuous, then the strains are continuous. Additional 
constraints can then be used to impose inter -element strain continuity but 
these conditions are not required for admissibility. The solutions obtained 
with and without these additional constraints will often have the same mean 
error (Reference 3) but the solution with strain continuity will require solving 
fewer equations. When dealing with higher-order elements the constrained 
stiffness matrix may be less than one-half the dimension of the original. 


When higher-order discrete elements with distinctly different strain- 
displacement equations must be assembled, a clear understanding of the 
admissibility conditions is essential. Perhaps the best way to describe the 
possible pitfalls is with an illustrative example. Consider adjacent flat plate 
and cylindrical panel elements from the discrete element model of the pear- 
shaped cylinder shown in Figure 1. Using the Bogner, Fox, Schmit (BFS) 
plate and cylindrical panel elements (Reference 4), there are 1Z degrees of 
freedom per grid point per element, the three displacement components 
relative to a local curvilinear coordinate frame u 0 ', u^ and the nine gradients 
of these displacement components u°> j, n a , 2 , u a , 12 , u 3 , j, u 3 , 2 , u^, 12. 
There is a tendency to erroneously assume the plate displacement gradient 
components are equal the panel displacement gradient components at a com- 
mon node, in particular (u^, 2 )fp “ ( u ^> Z)cP- ^ substantial error in the 
•normal displacement component u^ is caused by this assumption as shown in 
Figure 2 for the pear-shaped cylinder under a uniform axial load. The 
admissibility conditions merely require (u a ) F p = (u°0cp and (u^,^ )pp 
= (u^, Q/)cp where the curvilinear coordinates have been parameterized such 
that da = dS in each element. This allows discontinuities in u^, 2 which of 
course must exist if we are to have continuous midsurface strains, 
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( 1 ) 
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It should be noted that it is not necessary to use Equation (I) as a constraint. 

The minimum potential energy theorem ensures that the Ritz procedure for 

2 2 

admissible , displacement fields will find (u , 2 ^-pP an< ^ • 2 ^CP suc ^ ^at 
equilibrium is satisfied in the limit. In this problem equilibrium implies 
strain continuity and Equation (1) can be used as in Reference 5 as an additional 
constraint which reduces the number of equations to be solved. If we modify 
the pear-shaped cylinder such that the flat panel material is different from the 
cylindrical panel material then Equation ( 1 ) cannot be used since equilibrium 
now requires a strain discontinuity. 


There are at least two other situations in the assembly of higher -order 
discrete elements that deserve attention. Consider the stiffened cylindrical 
panel shown in Figure 3 again modeled using the BF$ discrete elements. In 
this case, even though adjacent cylindrical panel elements have the same strain 
displacement equations and are made of the same material, there can be dis- 
continuities in u\ ^ caused by load transfer between the stiffener and the 
cylindrical panel. Consider next an error that can occur when constraining 
the three rotations to be equal between adjacent higher-order plate elements. 

Let x^ = constant be the common edge between two BFS plate elements and 
recall that the elastic rotations are one -half the curl of the displacement vector; 
Using a Cartesian coordinate system, the rotations are 


1 3 

2 U ’ 2 


1 3 

" 2 u ’ 1 


e 3 = i (u 2 , , - a 1 . 2 ) (2) 

To ensure the same displacement along the common edge requires \i a , ? 

3 3 L 

and u , ^ to be continuous. If, in addition, we now require 0 to be continuous 

■ 2 

this will imply u , ^ is continuous. This in turn implies the shear strains are 

continuous which is erroneous if the plates are of different materials. Admis- 

1 2 

sibility requires only that 0 and 0 be continuous. 

These examples illustrate the pitfalls that can occur in the assembly of 
high-order discrete elements when constraints are used that exceed those 
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necessary for admissibility. Unfortunately it is not practical to use only 
admissibility constraints when strain continuity or other grid point con- 
straints are valid. These constraints not only reduce the number of 
equations, they usually do not change the structure of the stiffness matrix 
(if it was banded it will remain banded) and in most cases they do not increase 
the mean error. The design of a new structural matrix assembly module for 
use with higher-order elements in N AST RAN must take these factors into 
account. 


STRUCTURAL MATRIX ASSEMBLY MODULE 

The structural matrix assembly module in Level 15. 1 of NASTRAN can- 
not process discrete elements with more than six degrees of freedom per 
grid point. A new module is required for higher-order elements that accounts 
for accurate and efficient design requirements. Assembly based on simply 
the admissibility conditions must be available as a default and grid point con- 
straints must be available for efficiency. The module should be able to 
assemble the existing general elements in NASTRAN with higher-order 
elements of different types. This suggests two categories of grid point 
degrees of freedom for each element; those directly involved in admissibility 
conditions, UA, and all others, UH. All grid point degrees of freedom (in 
element coordinates) for all general elements now in NASTRAN fall in the 
first category. The new grid point degrees of freedom, UH, are somewhat 
like scalar point variables except they are elastically coupled to all the other 
grid point degrees of freedom for an element. As a default value, the number 
of UH at a grid point is equal the sum of the UH associated with that grid point 
from each element connected to that grid point. This corresponds to assembly 
based simply on admissibility. Next, it is necessary to provide for grid point 
constraints that are linear equations,, usually identities, among the UH and UA 
at a grid point. This is analogous to multipoint constraint equations with all 
the degrees of freedom occurring at the same grid point. As a practical matter 
a unique identification scheme for the UH will be needed. The UA of course 
already are identified uniquely by component numbers 1 to 6. Also, as a 
practical matter, an automated grid point constraint generator is needed; one 
that could set all UH components equal for elements of the same type at a grid 
point. To fix some of these ideas consider a BFS plate element (Reference 4), 
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a CKLO plate element (Reference 5) and a CQDPLT plate element 
(Reference 6) all modeling the behavior of a plate having one common grid 
point. At the element level the CQDPLiT element has five degrees of free- 
dom per grid point, the BFS element has twelve and the CKLO element has 
twelve. These are listed in Table 1, divided into UA and UH degrees of 
freedom. 


Table 1. Grid Point Degrees of Freedom 


Element 

UAj 

UA 2 

ua 3 

UA. 

4 

ua 5 ua 6 

BFS 

1 

u 

2 

u 

3 

u 

3 

u ’ 2 

3 

-a, j 

CKLO 

1 

u 

2 

u 

3 

u 

3 

u * 2 

3 

-u. j 

CQDPLT 

I 

u 

2 

u 

3 

u 

3 

u ' 2 

3 

-u, 1 


UHj 

UH 2 

UH 3 

UH. 

4 

UH C 

UH 6 

UH ? 

1 

1 

1 

2 

2 

2 

3 

U>1 

u > 2 

u > 12 

a, 1 

u * 2 

u ' 12 

u ’ 12 

1 

1 

2 

2 

Q 

3 

3 

U ’ 1 

u > 2 

u > 1 

a, 2 

u 'n 

u ' 22 

U ’ 12 


3 

The two higher-order elements have elastic rotations about the x axes (c.f. 
Equation 2) but these are not UA degrees of freedom as described earlier. 

If the UA^ degree of freedom is removed with an SPC1 card there are 19 
degrees of freedom at the common grid point. If the elements are all of the 
same material then grid point constraint equations can reduce this to 14. 
Suppose there is a rod normal to the plate that transmits torsion. This does 
not add a degree of freedom since UA^ in this case is related to Ut^ and UH^ 
of the BFS element by Equation (2) or equivalently UH^ and UH.^ of the CKLO 
element. As is obvious from this example, admissibility cannot be determined 
for all combinations of general elements without the analyst making judgments 
about load paths in his structure. These decisions will be input via grid point 
constraint equations. An analogous situation now exists with frame structures 
when the analyst uses pin flags (cuts and releases) on the CBAR card to input 
his decisions about joints. 


PARAMETRIC DISCRETE ELEMENT GEOMETRY 

The initial geometric representation of a complicated structure is a 
formidable design problem but it is one that has been solved by the time a 
discrete element analysis is required. Some form of a geometric model 
(loft lines, offsets, etc. ) has been prepared and serves as a data base for 
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the anhlyst. Increasingly these models are computer generated and in 

several industries piecewise polynomial surface representation is now used 

(Reference 7). This form of surface repre sentation is the same as that used 

for parametric discrete element models and consists of patches that map two 

12 3 

parameters (£, q) into spatial coordinates (x (£, q), ,x (^, q), x (g, q ), on 

the midsurface of the discrete element. The patches are constructed such 
that the edges of the element coincide with constant values of the patch para- 
meters (£, rj) as Figure 4 illustrates. The data required to define a patch 
with curved edges is obviously more than the grid point coordinates of the 
corners. If bicubic Hermite polynomials (Coons 1 surface patches) are used 
then x L , x l , g , x 1 ,^ and x 1 , ^ are required at each corner. To a large 
extent the increased data per grid point is offset by a reduction in the number 
of grid points required to model the geometry but this is a separate issue. 

The immediate problem is how to efficiently introduce into NASTRAN the 
geometric data required by parametric discrete elements. This^data can be 
input as property data for each element or as a separate entity like grid point 
coordinates that can then be referenced on a broader basis by all elements. 

The first approach would require a minimum change to NASTRAN but could 
be very inefficient in that the same boundary data might be input over and 
over again, once for each element sharing a common edge. Primarily for 
this reason only the second approach will be considered further. 

The patches used for parametric discrete elements are almost always 
bivariate polynomials although other interpolatory functions are possible. 

These polynomials can be uniquely determined in several different ways, each 
related to the other by a linear transformation. A bicubic Hermite poly- 
nomial interpolate determined from corner coordinates and derivatives also 
can be uniquely determined by the coordinates of sixteen points (Reference 1) 
where four of these must be interior points for the x\ ^ „ Boundary curves 
can also be used to define patches. Mallet provides an excellent example 
(Reference 8) using piecewise cubic interpolation of a grid line to obtain 
continuity. This approach uses grid line data in much the same way grid 
point data is now used in NASTRAN. Another useful representation is the 
super patch that defines several patches within its boundaries. The super 
patch is constructed using spline constraints (Reference 9) and has been used by 
Timmer (Reference 10) to form bicubic patches for aerodynamic surface model- 
ing. Both the grid line modeling and super patch modeling offer the additional 
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benefit of continuity along the entire common boundary between adjacent 

patches. Patches derived with the spline constraints of Reference 9 also have 
Z 

C continuity and require far less input data. This suggests a simple way of 

constructing grid line data from a near minimal data base. The parametric 

slopes at the two end points, x 1 , ^ ^ and x 1 , ^ and the grid point identification 

numbers Gl y GZ, . . . . , GN of points on the line are all that's required to define 

the line. If cross derivative data is desired for Coons surface patches then 

x 1 , ■ and x* T are also required but only at the four corners of a super 
g M 1 9 * i IN 

patch. A prototype data card for generating a grid line from spline con- 
straints is shown in Figure 5. There are of course situations where spline 
constraints produce a wavy line that does not model the initial geometry well. 
In this instance the parametric slopes should be input for every grid point on 
the line. Although no mention of grid line parameterization has been made, 
it will be necessary to adopt some standard such as 0 5£sl between adjacent 
grid points. 
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BULK DATA DECK 


INPUT DATA CARD SPLIN1 GRID LINE SPLINE CONSTRAINTS 

DESCRIPTION : DEFINE END SLOPES AND INTERMEDIATE GRID POINTS FOR 
A TYPE I SPLINE CONSTRAINT EQUATION 


S PL INI 

GLID 

CD 

XI Cl 

X2C1 

X3C1 

XI CN 

X2CN 

X3CN 

a be 

+bc 

G1 

G2 


• 

• 

GN 





FIELD 

CONTENTS 

GLID 

GRID LINE IDENTIFICATION NUMBER 

CD 

IDENTIFICATION NUMBER OF COORDINATE SYSTEM 
IN WHICH THE PARAMETRIC SLOPES ARE INPUT 

XI Cl. 

PARAMETRIC SLOPE X 1 , ^ AT GRID POINT 1 

X3CN 

PARAMETRIC SLOPE X 3 , ^ AT GRID POINT N 

GRID POINT IDENTIFICATION NUMBERS OF POINTS 
ON GRID LINE GLID IN SEQUENCE 

Gl. G2, .. GN 

Figure 5 

Prototype Bulk Data Card for Spline Constraints 
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